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ABSTRACT 

We calculate the two-body interaction of spherical hadronic bubbles im- 
mersed in a gluon plasma at temperatures above the phase transition. Mod- 
eling the bubbles with the MIT bag we find that the two body potential 
is repulsive for all bubble-bubble distances considered. This implies that a 
static configuration of spherical hadronic bubbles in a gluon background is 
consistent with the bag model of QCD. 
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1 Introduction 



Considerable theoretical work [|T], |^ and many computer simulations ^ 
of high temperature QCD have focused on the order of the phase transition 
between quark-gluon plasma and hadronic matter. Quenched simulations 
convincingly showed the transition is first order but weak i.e., with a rel- 
atively small latent heat. When dynamical quarks are incorporated, the 
transition weakens as smaller values of the quark mass are being used, and 
recent calculations suggest that in the limit rriq the transition becomes 
second order 0. 

Most investigations modeling this transition in the context of the Universe 
about 10-20 /isec after the big bang have assumed that the QCD transition 
proceeds by mechanisms similar to those encountered in first order phase 
transitions in condensed-matter systems, such as the condensation of water 
vapor. The common picture is of a homogeneous medium of QGP which cools 
adiabatically as the universe expands; as the temperature drops just below 
To, which is of the order of 100-200 MeV, the QGP becomes metastable, 
and hadronic bubbles nucleate by local fluctuations. These seeds expand 
slowly while exchanging heat and baryon number with the rest of the universe 
[g 1,0- 

Recently, however, an exact bag-model calculation suggested that the 
high temperature phase of QCD is not a simple homogeneous medium of 
quarks and gluons [0, |12]. It was claimed that at temperatures just above 
To a single bubble of finite radius of hadronic matter in the plasma has lower 
free energy than the plasma with no bubbles at all. This means that the 
high temperature phase of QCD may be composed of a plasma "soup" with 
floating bubbles filled mostly by an ideal gas of pions. However, the evidence 
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for this phenomenon from lattice Monte-Carlo calculations |T^, |14|, |T5|] is not 
clear and the interpretation of the results is not too convincing. Without 
addressing to the issue of how well the bag-model describes QCD physics at 
high temperatures, we try to understand whether within the arena of the 
bag-model such mixed-phase structure is indeed possible. 



The bag model calculation ||Tl| was performed for a single spherical hadronic 
bubble immersed in a homogeneous medium of plasma. The first question 
which arose was whether this spherical configuration is indeed a minimum of 
the free-energy functional and not a saddle point. In particular, is a spherical 
bubble stable against small deformations of shape? Using rather straightfor- 
ward analysis it was shown that perturbations of its boundary indeed tend 
to grow, suggesting that a spherical bubble may be unstable [0 . 

The aim of this work is twofold: to see if bubble-bubble interactions may, 
on the one hand, support a finite density of bubbles, and, on the other hand, 
stabilize the walls of the bubbles, such that the spherical bubbles discussed 
TT| , |TB|| are indeed stable objects within the scope of MIT bag model. 
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The next section will be devoted to the explanation of the strategy of 
obtaining an estimate of the bubble-bubble interaction. In section |^ we shall 
show how the multiple reflection expansion may be applied to calculating 
the density of energy levels for the simplest two-bubble geometry, and some 
details of the calculation will be given in section ^. In section ^ the results 
of the calculation will be displayed and in section | we shall discuss the 
implications of the results and the validity of the calculation. 
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2 The Method 



We assume that the high temperature phase of QCD is composed of identi- 
cal spherical domains (bubbles) where the fields arc manifested as colorless 
excitations (mainly pions); the bubbles are surrounded by freely propagating 
gauge fields. 

In the early universe, the rate of expansion is much slower than the time 
scales of QCD. This suggests that the field configuration is, at any time, in 
complete thermal equilibrium and, specifically, that there are no temperature 
gradients in the medium. Accordingly, we view it as a static configuration. 
We further assume that the fluctuations in bubble-bubble distances are small. 

Notice that in the absence of baryon number, the bubbles' radii and 
the mean distance between them completely specify the configuration. The 
structure we have in mind is rather similar to that of a type II superconductor; 
the hadronic domains form in our case a three dimensional lattice within 
the quark-gluon plasma, whereas in the type II superconducting phase the 
geometry is essentially two dimensional. 

Since bubbles are color singlet objects, their interaction should stem from 
a mechanism similar to the one producing the Casimir effect at T = 0. Such 
forces usually fall off very quickly as the relevant distance of the problem is 
increased. This leads us to assume that as long as the distance d between the 
bubbles is not too small, the interaction between bubbles is dominated by a 
two body interaction, i.e., that three body and higher order contributions are 
negligible in all configurations which are important to the partition function. 
Upon completion of the calculation, a consistency check for this assumption 
will be made. The smallness of the interaction enables us to express the free 
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energy of a two-bubble configuration as 

FbbiR, d) = 2Fb{R) + 6Fbb{R, d) , (1) 

where Fh{R, d) is the free energy of a single bubble of radius R in infinite 
space, and 6Fhb{R,d) is a small correction to be derived. For brevity of 
notation we omit the explicit T dependence. We first calculate the shift of 
energy level density Sp{k) for each gluon, from which the free energy shift 
SFbb{R,d) is derived via 

6Fbb{R, d) =ngT dk 5p{k) In (l - e'^'^) (2) 

where = 8 is the number of gluons. 

In the next section we shall briefly expound the multiple reflection ex- 
pansion, the tool with which 6p{k) is isolated and estimated. 

3 The Multiple Reflection Expansion 

The multiple reflection expansion is a method of calculating the distribution 
of eigenfrequencies for the wave equation in a connected finite domain of 
arbitrary shape. It was developed for scalar and vector fields |T^, |TB[ and 



later for the much more involved case of Dirac fields [T^. It uses the time 
independent Green function formalism for obtaining an expansion in 1/k 
where k is the wavenumber. The general idea is to extract the density of 
eigenfrequencies p{k) from the discontinuity of the propagator along the cut 
in the complex k plane, 

2 



Z f 

p-yik) = — Im / d\ G Arr' , k) 
TX I Jv 



(3) 



where V is the volume in which the fields propagate. 7 is an imaginary part 
added to k which measures the width of a Lorentzian smearing function. 
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chosen in order to turn the discrete eigenvalue distributions that pertains 
to finite domains into continuous ones. 7 must be chosen to be as small as 
possible so that k lies as close as possible to the real axis. 

G±{rr', k) is the physical propagator which satisfies the differential equa- 
tion together with the boundary conditions. In our model the fields satisfy 
the Helmholtz equation inside V , 

(V2 + A;2)l=0, divl=0, (4) 

and dual-superconductor boundary conditions on S, the boundary of the 
domain, 

h- E = hx B = 0. (5) 

As no confusion may arise, we shall henceforth omit the explicit k dependence 
of G. 

For deriving expressions for the energy level density it is technically con- 
venient to augment the transverse propagator with a longitudinal part which 
propagates longitudinal modes, 

G = Gx + G|| . (6) 

Using G for calculating P'y{k) via Eq. d), this additional piece yields an 
extra contribution which will be subtracted at the end. Proceeding along 
the lines of [|TB[, we define for each point a on the boundary of the domain 



the operators 



= Pn{a)£^-Pt{a), (7) 
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where Pn{ct) and Pt(«) are projection operators on the normal Ua and on 
the tangent plane to the surface at the point a respectively. Augmenting 
also the free transverse propagator with a longitudinal component, 



cxp(/A 


I'- -'■'!) 


An ' 


r — r' 





Go — Go_L + Goy — IGq , Go — — - — ■ — , (8) 



enables us to expand G in terms of Gq, yielding a multiple reflection expan- 
sion: 



G(rr') = Go(rr') +2 / cia„Go(m)DaDaGo(ar') 

s 

+2^ / daa dap Go(rQ;)DaDaGo(Qi/3)D^D/3Go(/3r') + • • • . (9) 



The first term gives for the density of states the usual contribution pro- 
portional to the (infinite) volume V ^ and is therefore irrelevant to 5p{k) . 
When computing the free energy of the system it is subtracted off. The sec- 
ond term, which corresponds to a single reflection, vanishes identically upon 
performing the trace and subtracting the contribution of the longitudinal 
modes. 

For the two-reflection term the integrand is 



Go{ra) — — - — (1 — sm B) 



dua drtD 



dGo{ra)OGo{<\ff) 



GoiPr') + 



dGQ{ra] 



Go{a(5) 



dGoiPr') 



dn 



sin' 9 



G'o(/3r')(2-sin2^), 



(10) 



where 9 denotes the angle between the normal vectors to the surface S at 
the points a and /5, namely cos^ = • n/3. Upon integrating over a, /3 and 
taking r = r' the triple products of propagators in the first and last terms 
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give identical results. Collecting all terms, subtracting the contribution of 
the longitudinal modes, and performing the trace, we finally get for the two- 
reflection shift in the density of states 



4 The Geometrical Setting 

In order to calculate the two-bubble interaction alone, it is sufficient to con- 
sider a simple geometry of two spherical bubbles in infinite space. This is in 
practice similar to confining the two bubbles in a cavity, and taking the walls 
of the cavity to infinity. The actual application of this limiting procedure is 
considerably simplified by the fact that only the bubble-bubble interaction 
is of interest and not the artificial interactions of the bubbles with the walls 
of the cavity. Therefore, we do not have to take into account refiections from 
the walls, i.e., the integrations over the points a and /3 are performed only 
over the two spherical surfaces of the bubbles. 

Figure |l]a shows schematically the relevant contribution to the propagator 
G, where the limit r —>■ r' has already been taken. The points a and (3 are 
on different spheres and r is in y. 

For the free energy shift 6Fi,h of the two-bubble configuration relative to 
the free energy of two completely separate bubbles, one has to subtract an 
excluded-volume contribution. Had we computed the free energy of a single 
bubble in infinite space by using the multiple refiection expansion, the two- 
refiection term would be as depicted in Fig. |l|b: The points a and (3 lie 
on the surface of the bubble, and the coordinate r runs over all space, the 




- sin^ e 



dua drip 



(Go(m)Go(a/9)Go(/?r)) 



(11) 
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Figure 1: Schematic representation of the two kinds of contributions to the 
two reflection term. See text. 

interior of the single bubble excluded. When introducing another bubble, its 
interior is no longer included in the region of integration, and therefore the 
excluded-volume terms must be subtracted off. 



5 Results 

We have computed the shift in the density of states 5p-^{k) as a function 
of the common radius R of the bubbles and the distance d between their 
centers. This enables us to estimate the free energy shift associated with the 
interaction of bubbles by Eq. (0). 



Since for the equilibrium bubble radius (see Refs. |TT|, |T^) RT ~ 1 , 
it was sufficient to compute Sp-y{k) in the range 0.2 < kR < 10. Lower or 
higher values of k are expected to yield too small contributions in Eq. (|^). 
The convergence parameter 7 was chosen as small as the numerics allowed. 
The results which follow correspond to the choice of 7 = OA/R. Calculations 
for lower values of 7 require much more iterations to achieve similar accuracy, 
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but lead substantially to the same results. Physically, such a high value of 7 
may be viewed as effectively introducing a gluon mass into the calculation. 

Integrating Eq. ([TT]) we encounter seven- dimensional integrals. Exploit- 
ing the azimuthal symmetry reduces them to six-dimensional integrals. We 
were, therefore, forced to resort to numerical routines which sample the huge 



space. We used the computer code VEGAS [^, which we found to be very 
efficient for the task.. We ran the numerical program on a Cyber 920 Silicon- 
Graphics machine. For each value of d the program ran for about 30 hours, 
sampling about 2 • 10® points for each value of k. 

The integration was performed using the bispherical coordinate system 
(see e.g. [^) in terms of which S is trivially expressed as a union of the 
surfaces /i = — /io and /i = /io where /i is the radial coordinate; /io fixes the 
ratio between the bubbles' radii and separation via d/2R = cosh/iQ. The 
volume integration is then performed over the range —fiQ < fi < fio for (a) 
type terms, and Hq < n for (b) type terms (see Fig. |ip. 

In general, the limit r —>■ r' must be taken with caution, since for a simply 
connected boundary all four points a, f3, r and r' may become arbitrarily 
close. In our case, however, the boundary is not simply connected and for 
both (a) and (b) type terms, the four points may not coincide. Therefore, 
the superficial singularities of the integrand of (|lID , a r or (3 ^ r for (a) 
type terms, and a — /? for (b) type terms, are not actual singularities due 
to the vanishing of the measure. From the numerical point of view, it is then 
sufficient to avoid sampling the space of integration too close to subspaces 



where the value of the integrand of (11) becomes unbounded. 

Figure |^ is a typical plot of the shift in the density of levels as a function 
of k for R = 1 fm and d = 2.59 fm. As expected, for small values of k, 
5p{k) approaches zero. This feature is common to all combinations of R 
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Figure 2: Typical behavior of Sp{k) for fixed values of R and d. 



and d we explored. The numerical integration, however, fails to yield the 
correct decrease of Sp{k) for large values of k, where the integrand is highly 
oscillatory. The insufficiency of the sampling mesh for large values of k is 
also reflected in the large error bars, which are of the order of the estimate 
itself. This does not prevent us from obtaining a reliable estimate for 6F{k) 
since for large values of k the logarithmic factor in Eq. |^ is exponentially 
small. 

The free energy shift as a function of the bubbles' radius and their distance 
is shown in Fig. |^. The relative errors, which are not shown here, are typically 
around 10 to 30 percent, though for the largest values of d/ R where SF goes to 
zero the relative errors are higher. The most important feature is the increase 
of 6F for small values of d/R which is most pronounced for 0.5 fm<,i?^2 fm. 
The meaning of this is that bubbles tend to repel each other, i.e., bubbles 
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do not tend to coalesce. This result supports the picture suggested in Refs. 
|n| , of a phase of quark-gluon plasma populated with hadronic bubbles 
arranged in a lattice-like structure. 

Another important result is that the largest values obtained for 5F are of 
the order of 100 MeV, and typical values are of the order of 10 MeV. This 
is to be compared with the free energy shift Fb due to a single bubble [see 
Eq. ([!)], which is typically about 500 MeV at the equilibrium radius. This 
means that the presence of neighboring bubbles should have little effect on 
the equilibrium radius. The mean separation between bubbles, however, is 
determined both by Ff, and 6Ff,f,: On the one hand, the tendency of Ff, is to 
populate the plasma with bubbles as densely as possible. On the other hand, 
5F pushes them away from each other as far as possible. The net result of 
these two mechanisms may be a stable lattice of bubbles. 

Three limiting regimes are worth mentioning. For d — 2R <^ R the results 
of the calculation are not shown, since for this case the distance between the 
surfaces of the bubbles is very small, and the propagator Go will decrease 
considerably only after many reflections. In other words, we expect that 
many reflections should be important in such situation and the two-reflection 
term badly approximates SF. For the limit R ^ with d fixed, the shift in 
the level density tends to zero due to the vanishing of the surface 5* in the 
integral (^. For d/R ^ oo the vanishing of 6 is clearly exhibited; this is to 
be expected due to the 1/R behavior of the propagator. 

We are particularly interested in the behavior of SF{R, d) as a function 
of d for fixed R. For each temperature we fix the value of R to the radius 
that minimizes the free energy for a single bubble geometry (see Ref. ||11|| ). 
The cuts through the free-energy surfaces yield the curves shown in Fig. ^. 
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Figure 3: SF for as a function of R and d/R for a temperature T — 155MeV. 
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Figure 4: 5F for plotted as a function of d, the distance between the bubbles 
for various temperatures. For each temperature R is fixed to the value which 
minimizes Fb{R,t). 



14 



6 Discussion 



We have estimated the shift of free energy related to the interaction between 
two spherical bubbles of hadronic phase immersed in a gluon plasma. In the 
relevant regime of parameters, a repulsion between the bubbles was found. If 
quarks play only a minor role in this phenomenon, as was found in the single 
bubble calculation in Ref. this finding supports the picture of the Swiss 
cheese instability ||12| of the early Universe. As an approximation we have 
used the first non-trivial contribution of the multiple reflection expansion. 
Since this is essentially an asymptotic expansion in it is not easy to 
predict a priori the validity of such truncation. 

At first glance, truncating all higher-order reflections has some similarity 
to the dilute gas approximation: For large d the free propagator Gq decreases 
strongly for every trip from one bubble to another. In other words, when the 
bubbles are dilute only terms which correspond to a minimal number of trips 
between bubbles are important. This may be misleading, since such terms 
arise at any higher order of the expansion. 

The correction we have computed is indeed considerably smaller than F^^, 
the free energy of a single bubble in infinite space. But an explicit compar- 
ison with terms stemming from larger number of reflections is completely 
impractical, because of technical problems in calculating integrals of more 
than nine dimensions. One may, however, be able to compute successive 
terms by using completely different ways of performing this calculation: For 



example, one might iterate the integral equation 11.29 of Ref. [jT8| on a grid, 
and use this as an input for equation 11.23. 

Our results are only qualitative in nature; The error bars are considerable, 
and we have only a limited window in k from which we may get reliable 
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numbers. Therefore it seems to us premature even to use them to predict 
the structure of the lattice of bubbles we suggest. 

The repulsion we found, is intuitively exerted by a positive outside pres- 
sure on the walls of the bubbles. However, by regarding the free energy as 
a function of the bubbles' radii for fixed d, it may be immediately noticed 
that it contains not only a pure pressure term, cubic in but other terms 
as well. These terms are crucial to the issue of the surface instability of 
the bubble, namely its tendency to grow fingerlike structures, as suggested 
in [^]. Again, much more precise calculation has to be performed in order 
to answer whether the two-bubble interaction may eliminate or reduce the 
surface instability. 
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